
#######################################################################
# Replication material for "The Social Dynamics of Collective Action: #
# Evidence from the Diffusion of the Swing Riots, 1830-31"            #
#                                                                     #
# Aidt, Leon and Satchell (2020)                                      #
#######################################################################

rm(list=ls())

library(spdep)
library(foreign)
library(stargazer)
library(sphet)
library(plm)
library(glmmML)
library(splm)

gc()

riots_panel <- read.dta("Swing_data_R.dta")
riots.p <- pdata.frame(riots_panel, index=c("input_fid","week_number"))
pdim(riots.p)
rm(riots_panel)


############################################
# Table 1, column 7: Fixed Effects Poisson #
############################################

res4<-glmmboot(totriots ~ L1riots_near+L1totriots+factor(week_number),
               family=poisson, data=riots.p, cluster=input_fid, boot=50)

summary.glmmboot
summary(res4)



#############################################################
# Interpret Fixed Effects results (column 1 of main table)  #
# This code is based on that in Mummolo and Peterson (2018) #
#############################################################

reg <- plm(totriots ~ L1riots_near + L1totriots + factor(week_number), data=riots.p, index=c("input_fid"), model="within")
summary(reg)
summary(reg)$coefficients[1:2,]

#Gives you the data (in rows) used for the regression
est<-reg$model
est["input_fid"] <- riots.p$input_fid
x<-"L1riots_near"
b<-reg$coefficients[x]


x.resid <- plm(L1riots_near ~ L1totriots + factor(week_number), data=riots.p, index=c("input_fid"), model="within")$residuals

reg2 <- lm(est$totriots ~ x.resid)
summary(reg2)$coefficients[1:2,]


y.resid <- plm(totriots ~ L1totriots + factor(week_number), data=riots.p, index=c("input_fid"), model="within")$residuals

reg3 <- lm(y.resid ~ x.resid)
summary(reg3)$coefficients[1:2,]

range(est[,x])
sdx <- sd(est[,x])

sdxresid <- sd(x.resid)
sdxresid

sdx/sdxresid

((sdxresid-sdx)/sdx) * 100


bounds <- list(NA)
ranges <- NA
units <- unique(riots.p$input_fid)
est$x.resid2 <- lm(est[,x] ~ est[,"factor(week_number)"] + est[,"L1totriots"])$residuals

for (i in 1:length(units)){
  bounds[[i]] <- range(est[est[,"input_fid"]==units[i],"x.resid2"])
  ranges[[i]] <- bounds[[i]][2] - bounds[[i]][1]
}

(length(ranges[ranges == 0])/length(ranges)) * 100


hist(ranges, col="grey", main="Within-unit ranges \nof Treatment",
     xlab = "Within-incumbent ranges", ylim=c(0,8000), xlim=c(0,20),
     axes=F)
axis(1, seq(0, 20, by=10))
axis(2,seq(0,10000, by=2000), las=2)
abline(v=mean(ranges), lty=2)
abline(v=quantile(ranges, probs=0.95), lty=2)
text(x=8, y=6000, labels="mean")
text(x=18, y=6500, labels = "95th\npercentile")
arrows(x=10, x1=mean(ranges), y0=5000, y1=5000, length=0.1)
arrows(x=20, x1=quantile(ranges, probs=0.95), y0=5000, y1=5000, length=0.1)

sdx
b*sdx

sdxresid
b*sdxresid

###########################################################
# The first two of these values are reported in the paper #
###########################################################
mean(ranges)
b*mean(ranges)

quantile(ranges,probs=0.95)
b*quantile(ranges,probs=0.95)

max(ranges)
b*max(ranges)
